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TRANSIENT EFFECT OF LUBRICANT ON 


ELASTOHYDRODYNAMIC FILM THICKNESS 


BY K. L. WANG and H. S. CHENG 

Department of Mechanical Engineering 
and Astronaut ical Sciences 
Northwestern University 
Evanston, Illinois 


SUMMARY 

The inlet solution of Elastohydrodynamic lubricated rolling contact problem 
was obtained considering lubricants with transient viscosity. The effect of 
the viscoelastic retardation time of lubricant on the center film thickness 
was investigated. 

1. The effect of transient viscosity in response to a sudden pressure 
was found to be insignificant in determining the film thickness in elastohydro- 
dynamic contacts. 

2. For the transient effects to become important in film thickness cal- 
culation, the retardation time would have to be at least three decades higher 
than those suggested by Harrison and Trachman in reference 9. 
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INTRODUCTION 


In lubrication of concentrated contacts such as rolling-element bearings, 
gears, and cams, it has been found by recent work on elastohydrodynamic (EHD) 
lubrication that the contacting surfaces are usually separated by a con- 
tinuous oil film. The level of this film thickness in elastohydrodynamic (EHD) 
contacts can be predicted by EHD Theories developed by Grubin (ref. 1), Dowson 
and Higginson (ref. 2), Archard and Cowking (ref. 3), Crook (ref. 4) and 
Cheng (ref. 5). Similar to the hydrodynamic theories in journal bearings, 
the minimum film thickness in EHD contacts was found not only to decrease with 
load and increase with speed and shear viscosity but also to be affected 
strongly by the pressure-viscosity dependence of the lubricant. In fact, it 
is because of this drastic increase in viscosity at high pressures, that 
contacting surfaces are separated by the hydrodynamic action of the lubricant. 

With regard to the accuracy of predicting the film thickness, the present 
EHD theories is only limited to moderately heavy loads and moderately high 
speeds. Recent work (refs. 6 and 8), have shown that there still exist large 
discrepancies between the isothermal EHD Theories and X-ray experiments for 
heavily loaded contacts. The inclusion of heating effects in the inlet of 
EHD contacts (ref. 7) accounts for some of the discrepancies, but the thermal 
theory does not predict a load dependence as strong as that measured by X-ray 
experiments . 

In searching for other possible reasons for this discrepancy, Bell and 
Kannel (ref. 8) suggested that the use of pressure-viscosity coefficients based 
on static measurements is invalid, because the increase in viscosity due to 
pressure rise in the high speed and heavily loaded cases may not behave in 
the same manner as measured in the static experiment. They developed a 
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Grubin-type inlet EHD theory assuming a short time-delay in the rise of viscosity 
with pressure. However, in their theory the selection of the time-delay con- 
stant is completely arbitrary, and what rheological mechanism governing the 
time-delay constant for a particular lubricant has not been studied. 

More recently, Harrison and Trachman (ref, 9) proposed a Transient 
pressure-viscosity model which enables one to predict the effective viscosity 
in the contact as a function of time. Using this theory, they have shown that 
the calculated effective viscosity as a function of rolling speed correlates 
very well with that measured by Johnson and Cameron (ref. 10) in the friction 
experiments. 

The object of this work is to incorporate Harrison and Trachman* s tran- 
sient pressure-viscosity model into the isothermal EHD Theory developed by 
Cheng (ref. 6), and to ascertain whether this transient pressure-viscosity 
effect will have a strong influence on the film forming capability in heavily 
loaded EHD contacts. 
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TRANSIENT VISCOSITY 


Doolittle's Empirical Relation 

Viscosity is a measure of fluid resistance to deformation and it depends on 
the state of fluid. Doolittle (ref. 11) adopted the idea that shear viscosity 
depends on the free volume of the fluid, which is defined as the free volume is 
the space when the liquid is expanded to a state from the state of absolute 
zero temperature. If is the specific volume of liquid at absolute zero 

temperature and v is the specific volume at normal state, then the relative 
free volume is defined as 

V - V 

( 1 ) 

o 

By performing a series of experiments, Doolittle found the following 
empirical relationship between viscosity and relative free volume. 

Tl^ = A Exp(B/f) (2) 

or 

In T1 = B/f + In A (3) 

s 

where A and B are material constants differed for each different liquid and B 
is usually very close to unity. This simple relationship will be used in later 
analysis to calculate the viscosity for a given state of free volume. 

Free Volume Viscosity and its Relation with Shear Viscosity 

The liquid structure can be interpreted by assuming that it is composed by 
a large number of crystal-like group of molecules. These groups of molecules 
undergo continuous breaking and reforming. Also, the atoms which should be in 
the neighborhood of some other atoms could be missing and thus produce a hole 
in that plate. The presence of holes adds an additional structural contribution 
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to the volume response of liquid when pressure or temperature is changed rapidly. 
If the pressure or temperature is suddenly changed, the liquid volume will under- 
go contraction or expansion and all molecules will rearrange themselves and pro- 
ducing more holes or filling up some holes. The latter process takes time to 
reach a new equilibrium state. By means of this structural relaxation process, 
the state of liquid after changes can be determined only when time scale is given. 

In order to describe this time-dependent behavior of liquid volume change, 
the following two simple models (Fig, 1) are used. 



Model A Model B 

Fig. 1 Models for compressional viscoelasticity 

Model A is a generalized Maxwell element with one relaxation time constant 
and model B is a special Kelvin element. Model A is convenient to correlate 
with experimental results and model B is good for later mathematical analysis. 

In model A, when a constant deformation is imposed, the stress p(t) 
follows 
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p(t) = [k^ + exp(-t/T)] Yq 


( 4 ) 


\ . 


where t = — is called relaxation time and in which T]^ is the volume viscosity, 

K 2 is the difference of instantaneous bulk modulus and equilibrium bulk 

modulus K - 
o 

By setting t = 0 in the time dependent modulus in equation (4) , one can 
easily get the instantaneous bulk modulus K = K + K« • When t = », this time 
dependent modulus becomes the steady bulk modulus as can be seen in equation. 
In model B, if a pressure p^ is imposed at time t = 0, the volume creep 
Y(t) can be written as 


i + • eXp(-t/T)]}p^ 


( 5 ) 


T is called retardation time, defined by T = — where TL is the free-volume 

Kf t 

viscosity and is the free volume bulk modulus. 

The instantaneous bulk compressibility ^ can be obtained by setting t = 0 

K. 

CO 

in time dependent bulk compressibility of equation (5), Also, the reciprocal 

of the steady bulk modulus is equal to ^ by simply inserting t = » in 

a> f 

equation (5) . 

A comparison of the modulus between two models yields 


= — 4. i- 

K K K. 
o 00 f 


( 6 ) 


= K + K. 
00 o 2 


( 7 ) 


Apply oscillatory bulk deformation and pressure to both models, one can 
get the complex bulk modulus as a function of frequency. 

For model A 


K = K + K„(iOU) = K + K„ — 

o 2^ ' o 2 1+ iu)T 


( 8 ) 
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For model B 


i = + i 

K K„ ^ K^(l + iuiT) 


relate equation (6) and (7) 


f _ 


relate equation (8) and (9) 


I » N 




( 9 ) 


( 10 ) 


( 11 ) 


Thus, there are two fixed equations (equation (10) and (11)) governing 
the relationships between the parameters of these two models. 

By measuring the propagation velocity and absorption coefficient of ultra- 
sonic waves propagated through liquid, Litovitz and Davis (ref. 12) obtained 
a method for calculating volume viscosity They found that volume viscosity 

is direct proportional to shear viscosity T]^, and it has the same temperature 
and pressure dependence as the shear viscosity. Since free volume viscosity 
is proportional to volume viscosity T] for a given state of liquid by adapting 

equation (11) where assuming the ratio — is known, it can be concluded that 

^2 

the free volume viscosity is proportional to shear viscosity T|^. 


Transient Response of Shear Viscosity to a Single Pressure Step 

A method originally derived by Kovac (ref. 13) for solving bulk creep 
behavior will be used here to calculate the transient shear viscosity of fluid 
after a finite imposed pressure step. Following his analysis, liquid having 
initial specific volume v^ will change to final equilibrium voliame v^ if there 
is enough time for change. With given value of P, the governing equation by 
using model B is 
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Vi - Vz = v^P/K^ (12) 

If the instantaneous volume change is which is equal to v^P/K^, equa- 

tion (12) can be written as 


(Vi - v^) + (v^ - v^) = v^P/K^ + v^P/K^ 


(13) 


it follows 


V. 


V, 


2 ^ 
K. 


(14) 


VI 

The time dependent part of volume change in model B can be solved from 
the differential equation considering force balance in parallel spring and 
dashpot combination 


\ A V. - V 

P = — ;^ + K 

dt f 


(15) 


substitute the value of P in equation (14) into (15) 




dv 


dt 


(16) 


for a finite change of pressure, can*t be considered as a constant since 
11^ is a function of dependent variable v- The governing equation becomes non- 
linear and it is difficult to solve. However, it was assumed in a previous 
section that free volume viscosity Tl^ is proportional to shear viscosity 
and both depend on free volume in the Doolittle’s empirical equation 


In 71^ = In A* + B/f 


(17) 


where constant B remains the same and close to unity. Define a parameter s 
such that 


s = In 


(y) = 


(18) 
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where f« is the final relative free volume for imposed pressure P and TL is 
the final equilibrium free volume viscosity. Equation (17) can be written in 
terms of parameter s. 


exp(-s) ds = - — 

s(l - s £ 2 / 8 ) T 2 


(19) 


where T 2 is a retardation time defined as 



( 20 ) 


and it will be evaluated at final equilibrium state. The term sf 2 /B in equation 
(19) is much less than unity so that (1 - sf 2 /B) ^ can be expanded and equation 
takes the form 


ds + exp(-s) ds = ~ 
s ^ B To 


( 21 ) 


For a given value of P, this equation can be solved numerically for s and by 
the relationship 


thus 


exp(-s) 


exp(-s) (22) 

It will give the time-dependent transient shear viscosity for liquid subject 
to a single pressure step. 


Transient Response of Shear Viscosity to a Continuous Pressure Chanpe 

In EHD problems, the lubricant moving through the gap between two rollers 
will experience a continuous pressure change from atmospheric pressure up to 
4 X 10^ psi within a very short time. Since this externally applied pressure 
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is a continuous one instead of a instantaneous pressure jump, the analysis used 
in the previous section cannot be used here directly. However, by approximating 
the continuous pressure input as a series of pressure steps as shown in Fig, 2, 
the previous method for solving the transient shear viscosity can be used re- 
peatedly and successively within each single step. In this case, equation (21) 
can be written as 


exp(-s ) 

ds . + exp(-s.) 

J ^ 3 




(23) 


and 


\ = \ exp(-s ) 

®j2 J 


(24) 


where variables with subscript j means it belonging to the j th pressure step, 



Fig. 2 Approximating the Pressure Distribution by Pressure Steps 
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Sj for pressure step j varies from initial value s^^ to final value From 


relationship 






since viscosity is continuous between each adjointing steps 


\ \ 

s . . ^ ■ 1 £ 

jl J-l,f 


Equation (25) becomes 


= 1“ (t^) * 1" (t^) 


J-1,2 




= F(P., P. i) + s. , . 

J J-1 J-l.f 

thus, initial value of for jth step can be derived from equation (27) once 

s. ^ f found in the previous stage. Referring to Harrison and Trachman (ref. 9) 

retardation time for most oils can be expressed as a function of the equilibrium 

shear viscosity T| and the pressure as follows 

50 

1 (28) 

3.5 X 10 + 9P 

for the jth pressure step, it becomes 


'"j2 


50 Tl^ exp(of»P,) 


3.5 X 10 + 9P. 


finally, after substituting equation (27), (29) into equation (23) and approxi- 
mating dSj by 




one obtains , 
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( 31 ) 


(s 


jf 



T^ (3.5 X 10^ + 9Pj) 
50 exp (or P^) 


Equation (31) is solved for 


’jf 


by using Newton's method. 
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GOVERNING EQUATIONS FOR FILM THICKNESS 


In formulating the elastohydrodynamic equations, the following assumptions 
are used: 

1. The rollers, as shown in Fig. 3, are subject to pure rolling. 

2. The deformation is purely elastic. 

3. The Hertzian width is much smaller than the width of the disks and 
the side leakage is neglected. Also, the Hertzian width is small in 
comparison with the disk radius so that the deformation can be cal- 
culated by the half-plane solution. 

4. The lubricant is isothermal and the inertia of lubricant is negligible. 
Equations governing deformation and pressure are 


^ 2 *2 
h = h + 




r 

ttE' J 


In 


Ll 


1 1 -X 


P(l)d| 


(32) 




(33) 


In non-dimensional form, above two equations become 
u -Ti fH - 

■ \*2) ^ \ \ „3 J 


dp 

dx 


H 


(34) 


H = 1 + 


— 2 
16 P 

HZ 


I* 

\2 rr J 


-1 - 

» II - X 1 


“0 


(35) 


The dependence of the equilibrium viscosity on pressure is assumed to be 
of the Barus form 


71 = Tl exp(c^) (36) 

s 2 o 

it follows from Equation (24) that transient viscosity Tl^ becomes 
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(37) 


Tig = T|^ exp(oP - s) 

Density change as a function of pressure is assumed as follows: 

where is the ambient density, c and d are constants from ASME Report (ref, 14). 

Equations (34), (35), (36) and (38) coupled with equation (31) can be solved 
by numerical method outlined in Appendix B and C. 
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RESULTS AND DISCUSSION 


Typical numerical results were obtained for a run with a load parameter 
^HZ 0.012 and non-dimensional center film thickness equals to 10 

A typical value of G(G = 3000) is chosen to illustrate the transient effect 
of the lubricant. The resulting speed parameter U for this case is equal to 
1.0307 X 10 which is very close to the value obtained by Cheng (ref. 6) 
without considering the transient viscosity. The results of inlet film thick- 
ness and pressure distribution for this run are plotted in Fig. 4. The ratio 
of transient viscosity to equilibrium viscosity as a function of the inlet 
position is plotted in Fig. 5. As can be seen in this figure, the viscosity 
ratio remains very close to unity over most of the inlet region. This shows 
that for typical conditions encountered in an elastohydrod 3 m.amic contact the 
response of lubricant viscosity to pressure is almost immediate in the inlet 
region. Since the film formation of an elastohydrodynamic contact takes 
place almost entirely in the inlet region, the transient characteristics of 
viscosity produce little effect on film thickness. However, in the center 
region, where the pressure is high, the lubricant viscosity does not respond 
to the pressure rise immediately. Since the frictional force in an EHD 
contact is largely governed by the viscosity in the center region, the transient 
effects become significant in the EHD traction calculation, as shown by Harrison 
and Trachman (ref. 9) . 

In order to determine at what level of retardation time the lubricant 

viscosity effects will become significant, a set of arbitrary multiplication 
2 3 4 5 

factors M = 10 , 10 , 10 and 10 is introduced for T^. Results for load 

parameter from 0.003 to 0.012 and normalized center film thickness H 

nZ c 

from 10 ^ to 10 ^ are shown in Table 1 and also are plotted in Figs. 5(a) 

— 2 

to 6(d) as a function of the rolling speed U. It is found that for M = 10 
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Table 1. OBTAINED NUMERICAL DATA 


* * , 

H = h /R 


'■hz ■ 


Values of Multiplication Factor M 


10 


10 -^ 


10 


10 " 


5.02.2901x10"^^ 5.093904x10"^^ 

7.212976x10"^^ 7.455977x10"^^ 


0.003 


0.00001 

0.006 


0.012 


0.003 

0.000005 

0.006 


0.012 

1 

0.003 

* 0.000002 

0.006 

1 

0.012 


0.003 

0,000001 

0.006 


0.012 


1.035142x10"^^ 

1.101366x10"’-’ 

1.86o423xl0'^^ 

1.867464x10"’-^ 

2.717763x10"^^ 

2.755632x10"’^ 

3.956260x10"^^ 

4.078289x10"’^ 

4.991350x10"^^ 

4.983464x10"’^ 

7.454592x10"^^ 

7.487535x10"’^ 

1.097072x10"^" 

1.105953x10"’^ 

1.835414x10"^^ 

1.830417x10"’^ 

2.788145x10"’-^ 

8.794191x10"’^ 

4.038028x10"’-^ 

4.066247x10"’^ 


5 i 885855x10"’^ 

1.087195x13"” 

9^659420x10"’^ 

2.136174x10"” 

1.701148x10"” 

4.167230x10"” 

2.011805x10"’^ 

3.098025x10"’^ 

3.190901x10"’^ 

6.026983x10"’^ 

5.167988x10"’^ 

1.183219x10"” 

5.104673x10"’^ 

-13 

6.361610x10 

7.915690x10"’^ 

1.167518x10"’^ 

1.230753x10"’^ 

2.195556x10"^^ 

1.844844x10"’^ 

2.055017x10"’^ 

2.861793x10"’^ 

3.573406x10"’^ 

4.253334x10"’^ 

6.173326x10"’^ 







Figure 3 - Geometry of lubricated rollers. 
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Distance from Center of Contact, x 


Figure 4 TOPICAL INLET FILM THICKNESS AND PRESSURE DISTRIBUTION 
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FIGURE 5 VARIATION OF RELATIVE VISCOSITY ALONG THE ROLLER 
CONTACT 
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1 


1 



10 


-11 


Velocity, U 


(a) Multiplication Factor M = 10 


FIGURE 6 FILM THICKNESS AS A FUNCTION OF VELOCITY FOR VARYING 
VALUES OF CONTACT STRESS AND MULTIPLICATION FACTOR M 








Thickness 


KJmyJL^ 



G=3000 


10 


-13 


i 



Velocity, U 


(d) 


Multiplication Facter M=10 


5 


Figure 6 (cont'd) 
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the film thickness is not reduced significantly comparing to that without 
considering the effect of transient viscosity. Significant reductions occur 
as M increases beyond two decades. 

The ratio of center film thickness calculated with the transient effect 

to that without this effect are plotted as a function of multiplication factor 

M in Fig. 7. It can be seen that significant reduction of film thickness 

3 

begin to occur when the multiplication factor M approaches 10 . It is somewhat 
unlikely that the level of retardation time of the lubricants under typical 
EHD condition can reach values several decades higher than those predicted 
by Harrison and Trachman (ref. 9). 
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Film 

Thickness 


0.8 


H 


(H ) 

c s 


0.7 



0.6 - 

0.5 I _j l_ 

10 ^ 10 ^ 10 ^ 

Multiplication Factor M 

a) U = 10‘^, G = 3000 



FIGURE 7 RELATIVE FILM THICKNESS AS A FUNCTION_OF MULTIPLICATION FACTOR FOR 
VARYING VALUES OF CONTACT STRESS AND U 




Film 

Thickness 




SUMMARY OF RESULTS 


The inlet solution of Elastohydrod 3 mamic lubricated rolling contact problem 
was obtained considering lubricants with transient viscosity. The effect of 
the viscoelastic retardation time of lubricant on the center film thickness 
was investigated. 

1. The effect of transient viscosity in response to a sudden pressure 
was found to be insignificant in determining the film thickness in elastohydro- 
dynamic contacts, 

2, For the transient effects to become important in film thickness cal- 
culation, the retardation time would have to be at least three decades higher 
than those suggested by Harrison and Trachman in reference 9. 
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APPENDIX A 


a 

A 

A* 

b 

B 

c 

d 


NOMENCLATURE 

semi-major axis of an elliptical contact 
constant in Doolittle's relation 

constant used in relation for free voliime viscosity 
semi -minor axis of an elliptical contact 
constant in Doolittle's relation 
coefficient in density function 
coefficient in density function 

9 

16 

48 U/H*^ 


1/E' 




^j2 


G 

h 

h 

o 

h* 

h 

c 

(h ) 

h ? ® 

min 


H 


H 


^«c)s 


1 2 , 2 
1 1 - 1 - 

Young's Modulus for rollers 1 and 2 
fractional free volume 

equilibrium state fractional free volume 

equilibrium state fractional free voltime for pressure step j 
qE* 

film thickness 

inlet film thickness at x = -b 


reference film thickness at = 0, 

dx ’ 

center film thickness at x = 0 


* 

h = h 


center film thickness for the case without transient viscosity effects 
minimum film thickness 


h/h 


h*/R 


h /R 
c 




s/R 
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^,3 


grid point numbers for the x coordinate 


k grid point numbers at x = x 

a a 

n iteration number 

K complex bulk modulus 

used in Eq, (h) , (i) and (j) 

low frequency bulk modulus 
high frequency modulus 

bulk modulus associated with molecular rearrangement of free volume 
for jth pressure step 

j 

complex relaxational modulus 
K 2 high frequency value of 


P 



P 

q 


pressure 





Q see Eq. (g) 

R R^R^/CRj^ + R2) 

RjjR^ radius of roller 1 and 2 





final value of s. 

J 


in pressure step j 
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T. 

J 


time required for lubricant pass through jth divided region 


U 


\<“l * "2> 

2E’R 


Ui,U 2 velocity of rollers 1 and 2 
V specific volume 

v^ specific volume at zero absolute temperature 

v^ initial specific volume 

V 2 final equilibrium specific volume 


V. 

X 

X 

* 

X 

X 

a 

X 

or 

ot 

T1 

s 


instantaneous volume response 

coordinate along the film 

reference coordinate at 4^ = 0 

dx 

coordinate separating the inlet region into two subregions 
coordinate separating the outlet region into two subregions 
x/b 

coordinate at the termination of the film 
HZ 

pressure-viscosity coefficient 
shear viscosity of the lubricant 

free volume viscosity 

equilibrium state free volume viscosity 




j2 




equilibrium state free volume viscosity for jth pressure step 
shear viscosity for jth pressure step 
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^j2 




P 


P 

o 

* 

P 


P 


T 


T 


T 


2 


T 


J2 


equilibrium state shear viscosity 

equilibrium state shear viscosity for jth pressure step 


inlet viscosity 

volume viscosity 

density of the lubricant 
ambient density 

* 

density at x = x 


o/p„ 


relaxation time 


Ko 


retardation time ■— 




t2 


K. 


K 


f . 

2 


^1"'^2 Poisson's ratio of rollers 1 and 2 

% dummy variable for x 

Y see Eq. (f) 
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APPENDIX B 


NUMERICAL ANALYSIS 

The region interested is the inlet half of the contact zone, which can be 
further divided into two sub-regions as shown in Fig. 8. In the first sub- 
region, pressure distribution is obtained by direct integration of the Reynold's 
Equation with introduced dimensionless function q where 


q = 1 - z- 

equation (29) can be written as 

dq _ 48U (E ~ p*/~p 1 

dx'H*2 "P ^ 

it can be integrated 


q(x) 


48U 

= J 

H -» 


d(ln V ) 
s 

dP 



(a) 


(b) 


(c) 


For a given viscosity as a function of pressure, the pressure distribution can 
be obtained by solving the equation 


TIJP) = 


1 - q(x) 

In the second subregion, the pressure distribution can be obtained by 
solving the combined equations (29) and (30). 

.3 j_ r- /— 2 , f« f l~Z ~ I \ — 


(d) 


H" dp r 1 


T) dx 
s 


P(5) In 


d?) - . 0 

1-51 / - -1 


(e) 


where — 16 P^^ /H and = 48U/H In the discretized form, it becomes 


Y, = 0. 

k 
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Distance from Center of Contact 


FIGURE 8 DIVISION OF PRESSURE IN THE INLET REGION 
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Y, = 




- c. 


^ \ 1 ^ 
k-^ 


1 + c 


(^. 1) 

K 2 


1 2 


kf-2 


i ^ PjQ(k- ^J) - :^] 

1=1 3 5 Pi 1 


(f) 


'■^k- — 

K 2 


These are a set of n equations to be solved by Newton-Raphson Method for P . 

K 

Where Q(K- “^jj) are the quadrature formulae for the singular logarithmic Kernel 
(ref. 6). 

3 

Q(k- 5,j> - 1 1 [K„(k.J)+ KJK_,.J)- K^(K^ - l.j)] (g) 


m=l 


where 


V. 


K 


lOc.jj) - ^ <-3Vj - ^ - UjClnlujI - 1) 


(h) 


2 

K2(k3jj) (Vj +V.^2> +7^ 


(i) 


■SO'.ij) - IT <-''j - ^ 

^ j 


(J) 


also 


'j - 


"j ■ 

2 


u. 


Vj = -^ (ln|uj I - |) 
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u. u. + 2 

’j “ ■ "j +2 ^’1+2 ■ 


(k) 


along with above equations , a set of n equations based on n grid points between 
-00 <0 can be rewritten again 

[exp(-s^)/s^] ds^ + exp(-s^) ds^ = - 

and 


(1) 


= '^O - "j> 

The following are the outlines of numerical procedures for solving the 
governing equations : 

1. Given a set of H, G values 

ciZ, 

2. Assume a pressure profile for -oo (0 

3. Calculate H(x) for (x (0 

4. Calculate density, viscosity for -co (O 

5. Integrate the following integral in the first inlet region 


(m) 


i(x) - r (2^^) « 

for -CO <x (x 


H 


6. Calculate U 

- 

^ = 48“'T^ 

^ a 

7. Solve equations (f ) , (1) by Newton- Raphs on method. 

8. Check the convergence for pressure. If not, repeat calculating pro- 
cedure from step number 3. 

9. Final solutions are in the forms of U, P, and H. 
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APPENDIX C 


NUMERICAL PROGRAM 


The complete computer program coded in FORTRAN IV is listed in this 
Appendix for solving the Transient Viscosity EHD problem. 
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u 


PROURAM HANG ( I^PU T » OUT PU T » PUNCH , TAPE5=INPUT ♦ TAPE6=0UTPUT ) 
tHUOl 

NASA-tMU Inlet ULM for line contact for load up to 400»000 hsi 
dimension op l6;> ,A(3 o,3o) »C Oo) t'^UMA (6 o» *K. AR (2Q) tPAC (10) 

COMMON P (60) ,M (5j ) (60) *Pm^PA (20)#'5(6C?60) »U8A (20) iHSa (20) 

common VISD(Cu))tL£N{60) tOENU (60) *PLUH (20) *SA (60) *SMA (60> »OX (60) 
common VlSl (3b) » vIS2(35) .VIS3(35) ,VIS(6?) 
common AF A,Pm^:cJ *hSB*UB*EU*LN,NR*NWtXF»KU,KR 
COMMON Pl,P2,UP\/.WTA*ITfUBO 
C READ BASIC INPUT DATA 

NH=5 
N ^ = 6 

READ(NRf 1) 
wRITE(Nw,l) 

PEAO (NR»2)fgKUN 
OO 1000 NHR=1 »NRlN 

REAO(iJH»2) K(3» NA* KO* KF» KR* NKER*NAV1S 

REAO(NR.2) NS1» nS 2» NS3* NS4* NSb* NSbi NS7» NS8* NS9t NSIO 
R£AO(Nn(i2) ITH» ITP, ITE 
REA0(NH*3) EHSm* EFSP* EPSE 
KNF=kH -1 

H£AUl!'W*3) (QX (N) »K = 1 »KKF) 

A ( I ) =-5.0 
UO n = 2,kF 

9SI X IK) =A («-l ) ♦DA (A-1 ) 

HZAO (NH, J) (P (K) ,K=1 ,KA) 
v‘(RITE (Nw * A > IK» P(iO* K = l* KA) 

HKA=P (KA) 

UO 1(J6 i^ = KA«kF 
Ijb P (K) =S JHT ( 1 ,0-A (K ) *-»2) 

NKa=K A- 1 
TEMPsP (KA ) /PKA 
UO 250 K=ltKKA 
p (K) =TEMP*P (^) 

PI=3, 141593 

C HEAD LOAD* SPLcQ AND LU3. PARAMETERS 

HEAU(NR»3) (PlU 31N)» N=1, b) 

REaO(nr,3) (K Al ( K ) » N=1 » B) 

PEAU(NR,2) NriH, XPhZ'M 
HEaU(NP*3) (mSaIin), N = 1* NMm) 

HEAL)(NR*3) (PriirjA(N)» N=l, NPMZM) 

IF (NAVIS ,EO. ^) 00 TO lOl 

hEAU (NR, 3) (Vijl (K) »K = 1, 31) 
hRITE (N.f,24) 

write (Nw, 20) ( VlSl (!<) ,K=1 ,3i ) 

DO i>)2 K=l*31 
lo2 VIS2 (K) =ALOG ( VISl (K) ) 
l-ll COI-JTINUc. 

WRITE (Nrt,7) 

wRITE(nw,2) kg« kA, kO* KF» KR» ^IKER 
wRI 1'E (f'<w , 8 ) 

wRltE(Nw*2) ITm* ITP* ITE 
WRITE (i'JW,9) 

WRITE (Nw,?o) EPSH*EHSP, EPSE 
wRlT[-; (N.-i,i j) 

WRITE (Nw,2 0 > (a (K ) ,K = 1 *KF ) 

IrLNKtR .EO. w ) 00 TO 9l 
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c 


CALL KEHCAL 


NS I 

IFCNSl .£»J. 0) GC TO 92 

rtKITE 

wrtlTE<NWf20) (((jHK»J)« J=l» KF)» K=1» KO) 

GO TO 92 

91 kEAD<N 9*21’ < » J> »J=1 »Kf ) *K*1 tKO) 

92 w.UTE(Nw,12) 
v*9ltE (iNW,2o) IPuJOtN) ^NsltG) 

W91 TE (Nw, l3) 

wrti TE (N<'»20 ) (hSA i \ ) »N=1 *NHM) 

wPITE 

WPITE <Nw*20^ ‘PM/BA <i\) »N=1»NPHZM) 

00 lOOO KPHZ=lfNFH2M 
PEAl)(NH*2) (KAmP (•■•i) »N = 1 »NHM) 

FHZrispHZsA (i\kHZ) 

■•'.RI TE (N‘^« 1^) RHAd 

AFAsPL'JH ( 1 ) Or-riZd 

tOs^^LOH (2) <>RnZ.->«2.0/'’l 
rf;=PLUd (3) ‘*HM4d"2.'j/PI 
IF("JAVIS .EO. «) GO TO 93 
HTArPLUrf (A) 

Rl=rLUd (S) /KdZJ 
pE=P|_Ud (c) /Pi-.Zo 
opv/r (i>2-p 1 ) / J j, 

(jPv2 = op v“2 . J 

VIS3 ( 1) = { vIS'i (2) - ts2 ( 1) •►Ah A*f)PV) /0PV2 
VIS3(31) = (VIS2(31) tdTA»OPv -VIS.2 (30) ) /O^VE 
UO ln3 k=2«3j 

1-3 VlG3(f<;) = lVlS4;(iS'>l)-VTS2(K-l) )/OPV2 

NSl 

JF( .IS1 .F-2. <) ) 0^ TO 93 
i^R ^ TE ,25) 

WRITE (•■('(,20) (ViS3 (K) ,K=l ,31) 

WRITE (fL.), 3) nPv, L)RV2 

93 LO 1 NH=i, ,-jrt.v 

U O R R = 

EACT = eACC-‘) 
hSHshGa (f'i-i) 

KA = KA.AR <f\H ) 

f.RI TF. <NWi 1 a) Hod 
G = HlL|R ( I ) 

C 1 = 1 h e >.''*^''Zd<>'»4;/P 5w 
C3 = -+H,r-Zp5-i'»*2 
Ca=C1/RI 

U!i(,= (ri5rt»T . 75/ (1,26 *G*h»0.6»PhZH«'>(- 0,27) ) )<»»(l0./7,) 

R J T E ( , 2 3 ) (j d G 
lT-1 

l'>7 C'ALL hCAL(KO) 

C NSl 

IFINSl .EG. o) GC TO 109 
w R I T E ( w , 5 ) 

V.RI TE U» P(K)» K = l* KO) 

1U9 CALL UVD (1,kO,2,0) 

CALL VUIC <FaCT,KA) 

CALL DVD (1, KO»l> u) 
r5 = Jr.ri (KC) 
rvK>\rKA* 1 
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STOP 

1 F0KMAT(72H 

I ■ ) 

2 FO«MAT<l6l5) 

3 FOKMAF (8E 10,3) 

^ FOKMAT (7 (lX»l2, 1 a» z13,6) ) 

5 F0HMAT(/6h h(K) /) 

6 fOHmaT (/6h P(K) /) 

7 format ( 5h KG» Sri KA» Sri KOt 5H KF» 5H KR» SH NKER ) 

8 FORMATtSh ITHf 5 m ITP* 5tl ITE) 

9 format ( loH EPSh flori EPSP * IqH EPSE ) 

ll» FOKMAf(6h XM(K) ) 

11 FOKMAT(7H (J(K*J)) 

12 F0KMAT(8h PLUd(N)) 

13 FOHv)AT(7h riSA(N); 

14 F0PMAT*9H PH/!dA<N') 

15 FO«MAT(5h PMZd=»tl 3.6) 

16 FOHMAT (dh SUmA(k)) 

17 FOHMAT (4h SOA) 

Id FOHMAT<5 h mSo=»£13.6) 

Pv format ( 1X»9E13. 6) 

21 FORMAT (5el5, O ■ 

22 FORMAT!/ 5d IT=» I5t 5 m Ufl=,El3.6) 

23 FORMAT (/inri UdoHLL<IM=»El3.6/) 

2A FOR-lATt / 8rl Vi5i(K)/) 

25 fOk-1aT(/8h VlSj(rt) /) 

26 format (33R MUlTIP .ICATION FACTOR FOR TAU 2 = *E13.6) 

ENO 
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CALL OVD 1) 

TP=M (KA) -OS/DLN (KA) 

1F(TP .LT.O.O) GC TO 116 
00 115 K=1*KA 

115 h (K) =ii (^) -TP 

116 SUMA(l)='r.^ 

00 135 K=1»KA 

lFlK-1) 117*117,118 

117 Zl = «H(t<)-05/Li£f^(K) )/H(K)«*3 
GO TO 119 

118 21= (H (K) -OS/OLN (M )/H(K)<*«3 

119 lFlK-1) 133,133,131; 

i'-id SOMA (k)=SUMA(K-1)*0.^*(X(K)-X<K-1))*(Z1*Z2)«\/ISD(K) 

133 Z2=2l 
136 continue 

S0A=1 .0-1 . 0 /vIS (KA) 

OC) 1155 p\=l,KA 

S0 = S0A<»5U^'A (k) / SOMA (KA) 

SSS = Sm IK) 

IF(SO_i,0 ) 1(>5, 1^1,141 
141 wRITUNv^^aI liMi soma »N=1 

k'iRITF. (|.j,vlSO (N) *N=1»KA) 

.-/RITE (;4,LE,S (N) »N=1»KA) 

V.' R i r F. ( N w , 4 ) ( N , rl I N ) , N = 1 , K A ) 

'.vRlTF. (imw»20) os 
125 P (K) =PMU (Sa.SSj) 

CALL riCAL (KU) 

C NS3 

IFINS3 .EO, u> QC TO 126 

KK'»=K A-4 

vv R 1 T E ( Ki >^ » = ) 

viRTlF. (N'.i,A) (K»h<K) .K = KK‘+»KA) 

.vRlTE ('N.',,6) 

viWi I E (N'.'J » 4 ) (K» F(k)» K=1* KF) 

i'jRlTE (N<^, 16) 

WRITE (N-^fA) (^, SOMA (K) »K = 1 *K A) 

1R6 CONTINUt 

UO = SOA/ (C3<»SuMa (KA) ) 

IFl ir *GT* l) uc= tUU+UOP) *0»5 
N=KO-KA 

c NS6 

IF(NS6 .£0* O) OC TO 136 
■/.HITE (Nil, 3) ci, C3, C4, Ob, OS, SOA 
■«R1TF(NW,4) (K.,ViS(K) ,K=1,K0) 

wHiFE (n»s,4) (k, visn (K) ,K=1 ,K0) 

wRITe (N/(,4) (K,LEN (K) ,K=1 , KO) 

WHITE (iNrt»4) (KtOEND (K) »K=1 ♦ KO) 

136 COijTINUE 
KKA = KA+ 1 
kKU = KU- 1 
00 I 70 K=KKA,KU 
HH= (h (K) +H (K-1 ) ) «U.5 
K K = K - K A 

DO l6iJ J = KA, KkO 
JJ=J-kA+ 1 

IF(J.EO.KA) bO TC 158 

A(KK,JJ)= C3^'Jb<»^ .50 CA«» (Q (K» J) -fO (K-1 , J) ) 

PO To 1-58 
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I5a SPQ=0. 

DO 137 L=1»KA 

137 SPQ=SPU* <Q (K»U+(i<K-l»L))*PtL) 

A(KK*JJ>= C3«Uo*w.5*C^«SPQ/P(KA) 

I3ti IF(J.E'J.K) GO rO 156 
IFlJ.EU.K-l) GU 10 157 
GO To IbO 

156 S10N=1.0 
GO TO 159 

157 SI«N=-1.0 

159 A (KK. JJ) =A (KK» JJ> *HH**3/(X I K ) -X ( K-1 ) ) /VIS ( K) 

1* <- (P (K) -P (K-1) ) *V1SD(K) *0.5*SI(3N) -C3*UB*0«5*DS*OEnO(K)/ 

2. DEN(K>««2 

\6u CONTINUE 

c (KK) =-HH**3/ (X(K)-X(K-l) )*{P(K)-P(K-1) ) / V I S ( K ) ♦C3*UB» ( HH-OS / 
IDEN (K) ) 

170 CONTINUE 

C NS4 

1F(NS4 .EQ. u) GC TO 174 
DO I7I KK=1«N 
WRITE U)w» 4 ) kk *C (KK) 

171 wRITE<N«*4) (JJ,AIKK»JJ) *JJ=1»N) 

J 74 CALL MATli\<V (A, N, C* 1» OtT) 

wPlTE ('^W,3) OET 

C NS 5 

IF (iNSS.EG.y ) GO 10 190 
wRJ lE (NW,‘0 IKK, C(KK), KK=1* N) 

1 9|j PkA = P ( K A ) 
r \/ = 1 * 

KK0=K0-1 
00 I5O KsKA.KKO 
kk=k-ka+i 
IJP (K) =C (KK ) 

IF I A^S (OF (K) ) -LP5P) I76,l7b,l75 

175 CO = c,:) 

176 P (K) =P (K) +0P (K) 
inU CONIINUE 

KKA=KA-1 

uo 1r1 6=1» KKa 

181 P (K) =P (K) »P (kA) / rKA 

IF(CV.EU.l.n) bO TO 210 
IF (IT.GT . IIP) ‘Gt TO 999 
IT=IT+1 

C NS7 

1F(N57 .EQ. 0) GC TO 200 
810 WRITE (Nw, 6) 

WRITE (NV<,4) (k, P(K)* K = 1* KO) 

WRITE ('''«*5) 

WRITE (NW.a) ( Mh(K)» Ksl* K0> 

WRITE (Nw, 16) 

wRITE(N.V, 4) (K*SUf'M(K)»K = l»KA) 

2uO write (N w, 22) IT.LU 
wRitE (Nw*26> FAC<pl 
IF (CV.EG.1.0) GC TO 999 
UBP-UB ■ 

GO TO IO7 
999 CONTINUE 

j.iu continue 
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SUBROUTINE HCAL »KK) 

COMMON P(60) ♦H(6C) fX(60) »PBiiBA ( 20) » Q ( 60 * bO) ♦UdA(20) »HSA(20) 

COMMON V ISO (BO) 'CEN (6o) ♦DENO(60) ♦PLUB(20, »SA<6o) ♦SMA(6o) ♦DX(6u) 

COMMON ViSl (j5}\vlS2(35) »VIS3(35) ,VlS(60j 

COMMON AFA»PHZiJ»f-Sa»UBfED»EN»NR»NW»KF».KO*KR 

COMMON Pl*P2»0PVtBTA»IT» U8U 

PI=3,lAl5q3 

CI=16.<»'^HZ8<»«2 /hSb 

00 10 K=1»KK 

h (K) =0. 

00 ^ J=1»nF 

1 h <K) =M (N) +P (0) *0 IK» J) 

H(K)=1.0+Cl»(y,5«X(N)»*2-H(K)/Pl) 

10 COIJTINUE 
RElUHrN 
E •'ll) 
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SUlii^OUTlNF. Kt^^^,AL 

COMMCU P(60)»Ml6l;)»X(60) tPHilBA(20) »Q (60*60) »uaA(20) »HSA(20) 

COMMON VlSD(bO) »l;EN (60) ♦U£ND(60) »PLUH(20) ♦ SA ( 6q ) t SM A { 6Q ) » DX ( 6y ) 

COMi-iON VISI (JS) ,VIS2(35) ,VIS3(35) ,viS( 60 ) 

common AF A*Pri/j»hS9»UB*£D*t.N*NR*NW*KF*K0»KR 

common pi f P2»UP^^*dTA» ITfUtJiJ 

DO 1 1=1* KF 

DO 1 J=l, KF 

1 Q(I*J)=0.0 

KKF=kF-2 
DO B K=l, KO 

G (Kt 1 ) =0 .0 

F5=X (K) 

DO 8 J = 1 ♦ KKF * 2 
U=X (J)-F5 
u2 = a ( J + <:) -f-5 
AU = AFJS (U) 

AU2= A8S(U2) 

IF ( A'J ) 6 si , 5u 

‘’0 AU = ALOO(Au) 

81 IF(AU2) 52* 6. 52 
52 A(J2 = ALaC3 ( AU2) 

& 0J=X(J+1)-X(j) 
p2 = J. L.'**L)j 
UO = IJOU 
Ij2g = ij2<nj2 

l-K = LlGO(Ab-l,5)*>0*5 

Fk2 = uA|j<» (aij2.1,5)'*Q,5 

FKM = lj'» (rcx-UG/bio) -U2tt ( F k2-u2q/ o . 0 ) 

G(K*J)= < 1-J.o*FK«FK2) /2.0-FKA/F2) (AO-1.0) + Q(K»J) 

«(K,J+1)= (2,>.'»(FK + FK2)+^r.0-A|-KR/F2)/UJ 

U(K*J + 2)=((-rX-3.o<*FK2) ^ 2 . 0-r K8/F 2 ) /G j ♦U2* ( AU2- 1 , 0 ) 

« COMTINUt 

DO 300 1^=1* «.0 
no 3oo J=l* 

3oo 0 , j ) =G , j) _U (XK, j) 

Hf. TUP fj 
ENU 
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SUBWOUTINE VOTU (FaCT»KA) 

DIMENSION FiJ<60) ♦TAU<60) 

COMMON P(6.0) »rtl6y) ,X(6o) »PhZRA(2o) »Q(6oi6o) ♦UBA(2o) tHSA(2b) 

COMMON VISO(6o) fCEN (6o> ♦DE nU ( 6o) ♦PLUfl (2o) »SA(6o) »SMA (6q ) fOX (6^ ) 

common VISl (35) , V1S2 (35) , VIS3 (35) ,VIS(60) 

common AFA,PriZD*rSR ,UQ,£D»tN,N'^»NW,KF»KO,KR 

COMMON P1»P2»OPV»dTAi IT,UB(i 

IF(lT.EO.l) Ub=UbG 

IF (it.GT.i) bO To 1 

siw=y . 00001 

S = SIN 

SM=AFA*P (1) 

1 ImAXs 5^y 

UENM=1.0 + EN<»o.-w ,5/Pm28/ (1.0*eD«0.015/PHZB) 

FPb= 0«0i)l 

do 5 n=UkO 

F2(N)= U;tNM-Ut.!. (N) ) /den (N) 

TAIJ (in) = 0.08-»UX (N)/US*PhZB**2« (PlUR(7)/Ph28 + 9*P (N) ) /£XP ( AFA*P (N) ) 
TAU(N)= TAU(i \|) /fact 
5 COfJTiNOE 

'•vPITt (Nw»A) ( TaC (N' ) * N=1 »K0) 

IF (lT»Gr.l) Go TC 12 
wPITE (^‘iv»l8u) 

GO 10 13 

12 S=Sa(1) 

SM=SmA ( 1 ) 

13 00 loo n=1»KU 
1 = 1 

10 ES-EXP(-S) 
tSM=EXP (-SM) 

T1=£S<»(F2(N)*1o VS) 

T2=ESH" (F-i (N) *i . w/SM) 

PSI=0.5-» ( Tl* 1 2) >(S-Sf.') ♦TAU (N) 

UPSI=u.3« ( (-Tl-£3/S**2) » (S-SM) ♦ (T1+T2) ) 

IF (PSI.GT. 0.1.) GO TO <^0 

US = -PSI/(;PSI 

9 IF (Arts IDS) -EPS) 11»11,13 

11 IF (ArtSiPSI) -EPS) 95,95»15 
15 IF (I-IMAX) 16»300»300 

lb 1=1+1 
S=S+DS 
GO TO 10 
20 S = S«(),1 
GO TO"10 

95 IK(N-KO) 96*9o*100 

96 s=S+OS 
SA('-'I)=S 

SM=S+AFA« (P (N+1 ) -P (N) ) 

SMA (N) =bM 
100 CONTINUE 
3«B CONTINUF. 

VIS(1)=1,0 
VISO (^ ) =AFA 
(JO 115 K = 2,kU 
IF (K-xA) JlOf- iC*3ll 
31o Pk = P('N) 

SiN = SA (K) 
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GO TO 312 

311 PK=‘(P <K) (K-1) ) *0.5 

SK® (SA (K) *SA (K-l ) ) *0.5 

312 VIS(K)=£XP(APA*P^-SK)‘ 

115 VI3L)(K)= APA-(.'5A(K)-SA(K-1) )/(P(K)-P(K-l) ) 

NS^® ■ 

IF(NSa.EG.O) oO TO 120 
A FOP-IAT (7 (lx, l2, lx, £13, 6)) 
l5o F0^^MAT (flh SA(X)) 

7 FOWmaT (1X,E1;j.6,8X,E13.6) 

155 FOA.-IAT (26H F2 TAU tIJ' 

171 POflMAT (IOh OIVEAGENT) 

IdO fOHMA r (3 x,1mN»2x,1hI, 5x,2h£S* 1 1 X , 3hESMi lOX » 2 hT 1 , lOX , 2 hT 2, lOx , 
l3HPSli lOX.^HUPbl » 10X.2HOS1 12X» IHS) 
lai format lix,213,8 (1X,E12.5) ) 

3j 0 *^RITE (r^w*17lJ 
12 u CONtl.-gUE 

return 

ENU 
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SUHROUTINE OVO (K1 *K2» IS.KHAF) 

COMMUN P(60) »H(60) ,X(60) tPHiBAiZO) »Q(60»60) fUBAtZO) *HSA(20) 
COMi^CN VISU(t>0) »L/En (6o) »UENO(60) »PLUB(20) »SA(6o) ,SMA(6o> fOXteu) 
COMMON VISU3ii) »VIS2(35) .VI53(35) ,VIS(60) 
common AFA»PhZ(J*PSP »UR»£0»tN»NH»NW*KFtKO»KH 
common p 1 ♦ P2 » dp V » OT a ♦ I T ♦ OHG 
IF(IS.EU.l) 2jj,2o5 
ECO H>k=P(K 2 ) 

IF (KrtAF.eO.l) PKs (P(K2) ♦P(K2-1) )«0.5 
D£N(K2)=1.0*EN*P«/ (1.0+ED*PK) 

DENO <k 2) =cN/ ( 1.0*PK*ED) »*2 
GO TO 2ib 

EoS UO 210 K=K1 ,k2 
PK=P (K) 

IFtKHAF .£0.1) Fk= (P(K> ♦P(K- l) ) »0.5 
D£N(k>= 1 .0*£N*Pl^/ ( 1 .0*ED*PK) 

Eli/ UENU^K)=E^/(E.c + ^'^*i:.D)«»2 
i£l5 COWTINUc 
hF TURN 
END 
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n on o n n non n n r. 


SUHt-iOUTIiNE MATiMV ( A » N i ti t M » UETtR ) 

A(3u,3^),^^(3o),IPIvo(3o),PIVOI(3o) 

c inversion with accOmRanyimg soluiiOn Of linear EOUATI' 

DETER =1.0 
UO 20 J=1»N 
2u IPI VO ( J) = j 
00 55^ I=1»N 

search for Pivot element 

AMAX=C.O 
00 105 J=1»N 

IK (IPIVC(J)-l) 6 vj,1o5.6o 
DO loo K=1»N 

IF <IPIVCIK) -I) Bo* loOf t>00 
Bw IF lABS (AMAX)-AdS (A(JiK))) 85,100*100 
Hi3 iROwsJ 
ICOLU =K 
AMAXsA ( J*K) 
iiiO CONTINUE 
lio CONTINUE 

1PIV0‘ICCLU)=IPIV0<IC0LU) ♦! 

INTERCHANGE KUrtS ]0 PUT PIVOT ELEMENT ON DIAGONAL 
IF <IR0W-IC0LU) 140* 260* 140 

1A(> ueter =-deter 

UO 200 L=1*N 

amAXsA ( irow*l) 

A ( IR0w,L> =A ( ICUlO*L) 

Zi, A ( IC0LU*L > =AMAX 
AMAX=H ( IROW) 

H ( IROW) =a ( ICOLU) 

H ( ICOLU) sAMAx 
pivot (I ) =A (ICOLU, ICOLU) 
deter =UETER*P1 VCT (1) 

DIVIDE PIVOT Row BY PIVOT ELEMENT 

A ( ICOLU* ICOLU) =1 .0 
00 350 L=1 *N 

35u A ( ICOLU* L) =A I ICOLU *L) /PIVOT (I) 

B ( ICOLU) =B ( ICOLU) /PIVOT (I ) 

REDUCE NCN-PIVOT ROWS 

380 DO 550 Ll=l *N 

IF(Ll-ICOLU) 400* 550* 400 
‘♦oo AMAX = A(L1*IC0LU) 

A(L1*1COlU) =o;o 
DO 450 L= I * N 

45^ A(Ll*L)=A (L1*U-A (ICOLU>L)*AMAX 
B(L1)=B(L1)-B(1CCLU)*AMAX 
55(i continue 
60’i> RETIJRN 
END 
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FUNCTION PMU(OU>SSS) 

COMMON H(60 ) ♦H( 6 c) ,X( 6o) »PH2HA(2o) ♦Q(6o<6o) »U8A<20) tHSAIZO) 
COMMON V1SU(6U) , DEN (60) »DE'^0(60) , PLUB (2'> ) » SA (60) »SMA(60) iOX( 60) 
COMMON VlSl <35* »VIS2<35) »VIS3<35) *715(60) 

COMMON AFA»PM2b*FSB*UB»E0*tN,NR*NW*KF»KO*KR 
COMMON Pl*P2*OPV»l3TA* IT*U9(i 
PMU=(-ALQGU,0-QCi) ♦SSS)/AFA 

Return 

ENO 
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